clear;clc;

load Data_CRSPMonthly N T permno yyyymm SzCat Excd Price Volat ESPR;

load PI_Continuous;
LambdaTS_Cont_Linear = NaN(N,T);
LambdaTS_Cont = NaN(N,T);
NobsTS = NaN(N,T);
Nobs0TS = NaN(N,T);
for t = 1:T
    idx = st_dataset.ym == yyyymm(t);
    st_t = st_dataset(idx,:);
    [~, ia, ib] = intersect(st_t.permno,permno);
    LambdaTS_Cont_Linear(ib,t) = st_t.lambda1(ia);
    LambdaTS_Cont(ib,t) = st_t.lambda2(ia);
    NobsTS(ib,t) = st_t.nobs1(ia);
    Nobs0TS(ib,t) = st_t.nzero(ia);
end
clear st_dataset t idx st_t ia ib;

NASD = Excd==3;
Micro = SzCat==3;
Small = SzCat==2;
Large = SzCat==1;
LnPrice = log(Price); LnPrice(~isfinite(LnPrice)) = NaN;
ESPR = ESPR(:,:,1);
clear SzCat Excd Price;

warning off;

WinsorLevel = 1/100;
NWLags = 12;

for loop = 1:2
    if loop==1
        Lambdas = LambdaTS_Cont;
    else
        Lambdas = LambdaTS_Cont_Linear;
    end

    Betas = NaN(T,18,5); AR2 = NaN(T,5);
    for t = 2:T
        Y = Lambdas(:,t);
        X = [Micro(:,t-1) Small(:,t-1) Large(:,t-1) ...
            NASD(:,t-1).*Micro(:,t-1) NASD(:,t-1).*Small(:,t-1) NASD(:,t-1).*Large(:,t-1) ...
            Volat(:,t-1).*Micro(:,t-1) Volat(:,t-1).*Small(:,t-1) Volat(:,t-1).*Large(:,t-1) ...
            -LnPrice(:,t-1).*Micro(:,t-1) -LnPrice(:,t-1).*Small(:,t-1) -LnPrice(:,t-1).*Large(:,t-1) ...
            ESPR(:,t-1).*Micro(:,t-1) ESPR(:,t-1).*Small(:,t-1) ESPR(:,t-1).*Large(:,t-1)];
        
        idx = all(isfinite(X(:,1:12)),2) & isfinite(Y);
        x = X(idx,:); y = Y(idx);
        y = winsorize(y,WinsorLevel);
        for ctr = 1:4
            [Betas(t,1:ctr*3,ctr),~,~,AR2(t,ctr)] = myols(x(:,1:ctr*3),y,2);
        end

        idx = all(isfinite(X),2) & isfinite(Y);
        x = X(idx,:); y = Y(idx);
        y = winsorize(y,WinsorLevel);
        for ctr = 5
            [Betas(t,1:ctr*3,ctr),~,~,AR2(t,ctr)] = myols(x(:,1:ctr*3),y,2);
        end

    end
    ToDisp = NaN(31,5);
    for ctr = 1:5
        n = 3*ctr;
        BB = Betas(:,1:n,ctr);
        for i = 1:n
            Y = BB(:,i);
            Y(~isfinite(Y)) = [];
            [B,~,Tnwest,~] = myols(ones(length(Y),1),Y,NWLags);
            ToDisp(2*i-1:2*i,ctr) = [100*B; Tnwest];
        end
    end
    ToDisp(end,:) = 100*mean(AR2(:,1:5),'omitnan');
    ToDisp
    disp(' ');
end
warning on;
